#include <cmath>
#include <iostream>
#include "vector"
using namespace std;
const double eps=1e-12;
class Function{
    public:
	    virtual double operator()(const double &x) const  = 0;
        virtual double diff ( double &x) {
            return ((*this)(x+eps) - (*this)(x-eps)) / (2*eps);
    }
};

class Polynomial : public vector <double>, public Function{
public :
    Polynomial (const int & n = -1, const double & x0 = 0, const double & x1 = 0) {
        (*this).resize(n+1);
        if (n >= 0)  (*this)[0] = x0;
        if (n >= 1)  (*this)[1] = x1;
    }
    int time() const  {
        return (*this).size() - 1;
    }
    virtual double operator ()(const double & x)const {
        int n = time();
        double p = 0;
        for (int i = n; i >= 0; i--){
            p = p * x + (*this)[i];
        }
        return p;
    }
    // 求导函数
    Polynomial  diff() {
        int n = time();
        Polynomial  p(n-1);
        for (int i = n; i > 0; i--)
            p[i-1] =  (*this)[i] * i;
        return p;
    }
    // 求导数值
    virtual double diff(const double & x)  {
        int n = time();
        double p = 0;
        for (int i = n; i > 0; i--)
            p = p * x + (*this)[i] * i;
        return p;
    }
    // #数乘的重载
    Polynomial  operator * (const double& x) const {
        int n = time();
        Polynomial  p(n);
        for (int i = 0; i <= n; i++) {
            p[i] = x * (*this)[i] ;
        }
        return p;
    }
    Polynomial  operator *= (const double & d)  {
        return *this = *this * d;
    }
    
    // 多项式加重载
    Polynomial operator + (const Polynomial & p1) const {
        int  m = min(time(), p1.time()), M = max(time(), p1.time());
        Polynomial  p2(M);
        for (int i = 0; i <= m; i++) {
            p2[i] =  (*this)[i] + p1[i];
        }    
        if (p1.time() <= time())
            for (int i = m+1; i <= M; i++) {
                p2[i] =  (*this)[i];
            }
        else
            for (int i = m+1; i <= M; i++){
                p2[i] = p1[i];
            }
        return p2;
    }
    Polynomial  operator += (const Polynomial & p)  {
        return *this = *this + p;
    }
    
    
    Polynomial  operator * (const Polynomial & p1) const {
        int n = time(), m = p1.time();
        Polynomial  p2(n+m);
        for (int i = 0; i <= n; i++)
            for (int j = 0; j <= m; ++ j){
                p2[i+j] +=  (*this)[i] * p1[j];
            }
        return p2;
    }
    Polynomial  operator *= (const Polynomial & p) {
        return *this = *this * p;
    }
    
};


Polynomial  operator * (const double & x, const Polynomial  & p) {
    return p * x;
}

ostream & operator << (ostream &out, const Polynomial  &p) {
    int n = p.time();
    if (n == -1) {cout << 0; return out;}
    for (int i = 0; i <= n; i++) {
        out << p[i];
        if (i == 1) out << "*x";
        else if (i >= 2) out << "*x**" << i;
        if (i < n) {
            if (p[i+1] >= 0) out << " +";
            else cout << " ";
        }
    }
    return out;
}

class Newton_Interpolation {
private:
    vector <double> x, y;
    vector <vector <double>> f;
public:
    Newton_Interpolation(const vector <double>& x, const vector <double>& y) : x(x), y(y) {
        int n = x.size() ;
        f.resize(n);
        for (int i = 0; i < n; i++) f[i].resize(n);
        for (int i = 0; i < n; i++) f[i][i] = y[i];
        for (int j = 1; j < n; ++ j) {
            for (int i = 0; i < n-j; i++)
                f[i][i+j] = (f[i+1][i+j] - f[i][i+j-1]) / (x[i+j] - x[i]);
        }
    }
    Polynomial  getPoly() const {
        int n = x.size() ;
        Polynomial  p(0, 0), pi(0, 1);
        for (int i = 0; i < n; i++) {
            p += f[0][i] * pi;
            pi *= Polynomial (1, -x[i], 1);
        }
        return p;
    }
};


class Hermite_Interpolation {
private:
    vector <double> x, y;
    vector <vector <double>> f;
public:
    Hermite_Interpolation(const vector <double>& x, const vector <double>& y) : x(x), y(y) {
        int n = x.size() ;
        f.resize(n);
        for (int i = 0; i < n; i++) f[i].resize(n);
        vector <int> m(n);
        for (int i = 0; i < n; i++) {
            if (!i || x[i] != x[i-1]) m[i] = i;
            else m[i] = m[i-1];
        }
        for (int i = 0; i < n; i++) f[i][i] = y[m[i]];
        vector <double> fac(n);
        fac[0] = 1;
        for (int i = 1; i < n; i++) fac[i] = fac[i-1] * i;
        for (int j = 1; j < n; ++ j)
            for (int i = 0; i < n-j; i++) {
                if (x[i+j] == x[i]) f[i][i+j] = y[m[i]+j] / fac[j];
                else f[i][i+j] = (f[i+1][i+j] - f[i][i+j-1]) / (x[i+j] - x[i]);
            }
    }
    
    Polynomial  getPoly() const {
        int n = x.size() ;
        Polynomial  p(0, 0), pi(0, 1);
        for (int i = 0; i < n; i++) {
            p += f[0][i] * pi;
            pi *= Polynomial (1, -x[i], 1);
        }
        return p;
    }
};

Newton_Interpolation Grid_Interpolation (const Function  & f, const double & l, const double & u, const int & n) {
    vector <double> x(n), y(n);
    for (int i = 0; i < n; i++) {
        double xi = l + (u - l) * i / (n-1);
        double yi = f(xi);
        x[i] = xi, y[i] = yi;
    }
    return Newton_Interpolation(x, y);
}

Newton_Interpolation Chebyshev_Interpolation (const Function  & f, const double & l, const double & u, const int & n) {
    vector <double> x(n), y(n);
    for (int i = 0; i < n; i++) {
        double xi = (u + l) / 2 + (u - l) / 2 * cos(M_PI*(2*i+1)/(2*n));
        double yi = f(xi);
        x[i] = xi, y[i] = yi;
    }
    return Newton_Interpolation(x, y);
}